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(Dated: February 2, 2008) 

The stress-strain relations and the yield behavior of a model glass (a 80:20 binary Lennard- Jones mixtures jjl) 
is studied by means of molecular dynamics simulations. In a previous paper |2] it was shown that, at tempera- 
tures below the glass transition temperature, Tg, the model exhibits shear banding under imposed shear. It was 
also suggested that this behavior is closely related to the existence of a (static) yield stress (under applied stress, 
the system does not flow until the stress a exceeds a threshold value Oy). A thorough analysis of the static 
yield stress is presented via simulations under imposed stress. Furthermore, using steady shear simulations, the 
effect of physical aging, shear rate and temperature on the stress-strain relation is investigated. In particular, we 
find that the stress at the yield point (the "peak"-value of the stress-strain curve) exhibits a logarithmic depen- 
dence both on the imposed shear rate and on the "age" of the system in qualitative agreement with experiments 
on amorphous polymers |3, 4] and on metallic glasses |5, 6]. In addition to the very observation of the yield 
stress which is an important feature seen in experiments on complex systems like pastes, dense colloidal suspen- 
sions 01 and foams [3], further links between our model and soft glassy materials are found. An example are 
hysteresis loops in the system response to a varying imposed stress. Finally, we measure the static yield stress 
for our model and study its dependence on temperature. We find that for temperatures far below the mode cou- 
pling critical temperature of the model (Tc = 0.435), dy decreases slowly upon heating followed by a stronger 
decrease as Tc is approached. We discuss the reliability of results on the static yield stress and give a criterion 
for its validity in terms of the time scales relevant to the problem. 

PACS numbers: 64.70.Pf,05.70.Ln,83.60.Df,83.60.Fg 
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I. INTRODUCTION 

Despite the large diversity of their microstructures, the so 
called soft glassy materials |9] like pastes, dense colloidal sus- 
pensions, granular systems and foams exhibit many common 
theological properties. Once in a glassy or "jammed" state, 
these systems do not flow, if a small shear stress is applied on 
them. For stresses slightly above a certain threshold value (the 
yield stress, o-y), however, they no longer resist to the imposed 
stress and a flow pattern is formed f^.wHUO. 

Let us illustrate this behavior using results of simulations to 
be described in more detail in later sections. Figure^shows 
CAnax, the maximum velocity in the system measured close to 
the left wall during simulations of a 80:20 binary Lennard- 
Jones (LJ) mixture 1 1 ] while applying a constant shear stress 
to the left wall (see sectionlUfor more details on the model). 
The applied stress is increased stepwise by an amount of da = 
0.02 every 4000 LJ time units and Umax is measured between 
two increments of the stress (note that, as seen from the inset 
of the same figure, this time is long enough in order to also 
determine the velocity profile, u{z), accurately). 

It is seen from Fig. ^ that, for stresses a < 0.6, C/max is 
hardly distinguishable from zero. In particular, it is much 
smaller than the thermal velocity of the wall, [/*ermai _ 
^/T/M w 0.0236 {M = 360 is the mass of the wafl and 
T = 0.2 the temperature). Thus, at these stresses, the sys- 
tem remains in the jammed state and resists to the drag force 
transmitted to it by the left wall. However, as the stress is fur- 
ther increased, a remarquable change in the system mobility 
is observed. The system starts to flow and Umax increases by 



more than two orders of magnitude. 

An inspection of the corresponding velocity profiles illus- 
trated in the inset of Fig.^reveals a further feature related to 
the yield stress, namely that, once the applied stress exceeds 
the yield value, the whole system fluidizes and the velocity 
profile is practically linear (velocity profiles corresponding to 
cr < 0.6 fluctuate around zero and are not shown in the inset). 

On the other hand, in experiments upon imposed shear rate, 
shear thinning is observed The apparent viscosity, 

defined as the average stress divided by the average overall 
shear rate, ?7app = a/jtoi, decreases with increasing jtot (in 
the case of a planar Couette-flow with wall velocity and sep- 
aration U^aii and L^, for example, 7tot = U^aii/L^). Fur- 
thermore, over some range of shear rates, the system sep- 
arates into reg ions with different velocity gradients (shear 
bands) ElHlIl- 

Whereas the shear thinning is commonly attributed to the 
acceleration of the intrinsic slow dynamics by the external 
flow (the new time scale, 1 / 7tot, is much short erjhan the typi- 
cal structural relaxation time of the system) |9, 16, 17, 18, 191], 
the origin of the shear bands still remains to be clarified. In 
some cases, this shear-banding phenomenon can be under- 
stood in terms of underlying structural changes in the fluid, 
analogous to a first order phase transition. Examples are 
systems of rod like particles, entangled polymers or surfac- 
tant micelles where the constituents (rods, polymer or surfac- 
tant molecules) gradually align with increasing shear rate thus 
leading to a coupling between the local stress and the spatial 
variation of the velocity gradient li2Qt i21i] . In the case of soft 
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applied stress 



FIG. 1: The maximum velocity in the system, C/max. measured in 
the layer of closest approach to the left wall during simulations of a 
binary Lennard- Jones glass [T = 0.2 (< Tc — 0.435)] at imposed 
stress. The stress is increased by an amount of da — 0.02 once in 
4000 LJ time units and UmnK is measured between two subsequent 
stress increments. The horizontal dotted line marks the thermal ve- 
locity of the wall. Note the sharp increase in i7max when changing 
the stress from a = 0.6 to 0.62. Inset: rescaled velocity profiles, 
u{z) /[/max, measured during the same simulations for stresses in the 
flow regime (for which [/max > [/*=™^'). Obviously, once the flow 
sets in, a linear velocity profile is formed across the system. 



range of shear stresses, a situation encountered in several com- 
plex fluids 1 2^ . This phenomenon should thus be generic for 
many soft glassy materials. 

In this paper we present an extensive study of the stress- 
strain relations and yielding properties of the present model. 
The report is organized as follows. After the introduction of 
the model in the next section, results on the system response 
to an imposed overall shear rate are presented, and the effects 
of physical aging, shear rate and temperature on the stress- 
strain curves are investigated. In section HV1 the response of 
the system to imposed stress is studied. The measurement of 
the static yield stress in the subject of section|V] A summary 
compiles our results. 



II. MODEL 

We performed molecular dynamics simulations of a generic 
glass forming system, consisting of a 80:20 binary mixture of 
Lennard-Jones particles (whose types we call A and B) at a 
total density of p = pa + Pb = 1-2. A and B particles interact 
via a Lennard-Jones potential, Uu{r) = '^iapli'^ap/'r'Y'^ — 
((Tc^/r)^], with q;,/3 = A,B. The parameters eaa, caa ™d 
toa define the units of energy, length and mass. The unit 
of time is then given by r = ctaa V^^^T^^a- Furthermore, 
we choose eab = I-Seaa, ebb = O.Seaa, ctab = O.Sctaa, 
<^BB = O.SSctaa and rris — rri/^. The potential was truncated 



glassy materials, however, no such changes are evident, and 
coexistence appears between a completel y steady region (zero 

shear rate) and a sheared, fluid region |8, "ullUlllii. 

It was shown in a previous work |2] that a model of 80:20 
binary Lennard-Jones glass Qj] also exhibits the shear band- 
ing phenomenon. Furthermore, a link was suggested between 
the occurence of shear bands and the existence of a static 
yield stress in the system. It was found that [see Fig. |3 the 
yield stress is larger than the steady state stress measured in 
a steady shear experiment in the limit of the zero shear rate, 
(Ty > lim.y,_^|^o cr. It was then suggested that, a shear-banding 
could be expected for shear rates, for which cr (7tot) < <Jy'- 
as the flow is imposed externally (by moving, say, the left 
wall) the formation of a flow pattern is unavoidable. On the 
other hand, it follows from cr(7tot) < fy that, some regions 
in the system are "rigid enough" to resist to the flow-induced 
stress whereas other regions undergo irreversible reaiTange- 
ments more easily |24]. Hence, whereas the details of the 
"nucleation" and growth of a heterogeneous flow pattern may 
depend on the initial heterogeneity in the "degree of jam- 
ming" |25], "free volume" |26] or "fluidity" |27, 28] at the 
beginning of the shear motion, its very origin lies in the pos- 
sibility of resisting to the shear-induced stress, i.e. in the ex- 
istence of a static yield stress. 

Therefore, although it does not solve the problem of the se- 
lection between the two bands, the existence of a static yield 
stress is at least consistent with the coexistence of a jammed 
region and a fluidized band: once the yield stress Cy is added 
to the flow curve, the shear rate becomes multivalued in a 




FIG. 2: The shear stress versus imposed shear rate, 7tot, under ho- 
mogeneous flow conditions at T = 0.2. The square on the horizontal 
axis marks the yield stress measured in imposed stress simulations of 
a planar Couette cell [see the dramatic change in [/max at cr = 0.6 in 
Fig.Q see also Fig. ll7l . Under imposed shear, if the corresponding 
steady state stress falls below the horizontal dotted line, a heteroge- 
neous flow can be expected, whereas in the opposite case the flow 
will be homogeneous. The vertical dashed line marks the shear rate 
on the boundary of these two flow regimes. Note that the yield stress 
shown here is a lower bound for Oy (see the solid line in Fig. I17t 
and thus is smaller than the value used in ^2*1. However, as a com- 
parison with Fig. 3 of Ref. |2] shows, the estimated 7tot-range for 
heterogeneous and homogeneous flow regimes is hardly changed by 
this modification. 
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at twice the minimum position of the LJ potential, = 2.245. 
Note that the density is kept constant at the value of 1.2 for 
all simulations whose results are reported here. This density 
is high enough so that the pressure in the system is positive 
at all studied temperatures. The present model system has 
been extensively studied in previous works dlisllallal and 
exhibits, in the bulk state, a computer glass transition (in the 
sense that the relaxation time becomes larger than typical sim- 
ulation times) at a temperature of ~ 0.435 1 1]. Since our 
aim is to study the interplay between the yield behavior and 
the possible flow heterogeneities, we do not impose a constant 
velocity gradient over the system as done in Ref. 1 18], where 
a homogeneous shear flow was imposed through the use of 
Lees-Edwards boundary conditions. Rather, we confine the 
system between two solid walls, which will be driven at con- 
stant velocity. By doing so, we mimic an experimental shear 
cell, without imposing a uniform velocity gradient. 

We first equilibrate a large simulation box with periodic 
boundary conditions in all directions, at T = 0.5. The system 
is then quenched to a temperature below Tc, where it falls out 
of equilibrium, in the sense that structural relaxation times are 
by orders of magnitude larger than the accessible simulation 
times. On the time scale of computer simulation, the system 
is in a glassy state, in which its properties slowly evolve with 
time towards the (unreachable) equilibrium values (aging, see 
Fig.|5ji. After a time of t = 4.104 [2.10'' MD steps], we create 
2 parallel solid boundaries by freezing all the particles outside 
two parallel xy-planes at positions z^aii — ^l^zl'^ (Lz — 40) 
[see Fig.EJ. For each computer experiment, 10 independent 
samples (each containing 4800 fluid particles) are prepared 
using this procedure. Note that the system is homogeneous in 
the a;?/-plane (L^ — Ly ~ 10). We thus compute local quanti- 
ties like the velocity profile, the temperature profile, etc. as an 
average over particles within thin layers parallel to the wall. 

The amorphous character of our model is clearly seen by an 
analysis of the packing structure, i.e. the radial pair distribu- 
tion function. Figure|3]shows the various kinds of radial pair 
distribution functions which can be defined for a binary mix- 
ture: gap is the probability (normalized to that of an ideal gas) 
of finding a particle of type a at a distance r of a /3-particle 
(q!,/3 £ {A,B}). In order to demonstrate that the system keeps 
its amorphous structure at temperatures far below the glass 
transition temperature of the model, we show the mentioned 
pair distribution functions at two characteristic temperatures, 
one in the supercooled state (T = 0.5 > Tc = 0.435) and one 
at T = 0.2. As seen from Fig. |3j the maxima of gap are more 
pronounced at lower T. However, no sign of crystallization or 
long range positional order is observed as the temperature is 
lowered through the glass transition. 

The mentioned insensitivity of the static structure to the 
glass transition must be contrasted to the fact that, at tempera- 
tures slightly above Tc, the system can be equilibrated within 
the time accessible to the simulation whereas this is no longer 
the case for temperatures significantly below Tc. At T = 0.5, 
for example, the time necessary for an equilibration of the sys- 
tem is of order of a few hundred Lennard-Jones time units (not 
shown). For T = 0.45, the equilibration time rises to a few 
thousands whereas at T = 0.2 the system is not equilibrated 
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FIG. 3: The radial pair distribution functions, g^p (a, /3 G { A,B}) 
at two characteristic temperatures of r = 0.5 (supercooled state) and 
T = 0.2 (glassy state, note that the mode coupling critical temper- 
ature of the system is Tc = 0.435). Note that curves belonging to 
T = . 5 and T = . 2 are qualitatively similar (no further peaks occur 
while cooling below T). Thus, the system maintains its liquid-like 
(amorphous) structure at temperatures far below T^. 




FIG. 4: A snapshot of the system at a total density of p = pA + Pb = 
1.2 (pA ~ 0.96 and pB = 0.24). The walls (darker particles) are 
made of the same types of particles as the fluid itself. They are dis- 
tinguished from the inner particles in that either they have no thermal 
motion or they are coupled to equilibrium lattice sites by harmonic 
springs thus preventing their diffusion. 



even after 2 x 10^ LJ time units. At this temperature, time 
translation invariance does not hold and the dynamical quan- 
tities depend on two times: the actual time, t, and the waiting 
time t„. Here, is the time elapsed after the temperature 
quench (from T = 0.5toT = 0.2) and the beginning of the 
measurement. 

This behavior is illustrated in Fig.|5j where the mean square 
displacement (MSD) of a tagged particle is shown at a temper- 
ature above Tc (T = 0.45) and at T 0.2 (far below Tc) for 
various waiting times. The figure nicely demonstrates the es- 
tablishing of the time translation invariance (TTI) at T = 0.45. 
Here, t„ = corresponds to a change of temperature from 
T = 0.5 to T = 0.45. As expected from the fact that T = 0.45 
belongs to the supercooled (liquid) state, with increasing wait- 
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ing time, the MSD converges towards the equiUbrium curve 
reaching it after about 4000 Lennard-Jones time units. It is 
worth noting that the waiting time at which the TTI is recov- 
ered roughly corresponds to the time needed for the MSD to 
reach the size of a particle. 

At T = 0.45, the equilibrium curve for the MSD exhibits 
the well-known two step relaxation characteristic of super- 
cooled liquid: for short times (t <C 1), free particle motion 
with thermal velocity is observed ([r{t + t^) — r(tw)]^ — 
[v^^tY — ik-^Tt^). The free (ballistic) motion ends up in a 
plateau thus indicating the (temporal) arrest of the tagged par- 
ticle in the cage formed by its neighbours. Akeady after a 
few hundred LJ time units, the plateau is gradually left and 
the MSD crosses over towards a linear dependence on time 
(diffusive regime). This is indicative of cooperative relaxation 
processes leading to the final release of the tagged particle 
from the cage (cage relaxation). 

At r = 0.2, however, the situation is completely different. 
Here, TTI is not reached on the simulation time scale. Even 
after a waiting time of 10^ LJ time units, the MSD continues 
slowing down without reaching a steady state. The slowing 
down of the dynamics with t„ also has a direct consequence 
on the life time of the cage. Figure |5] shows that, as i„ in- 
creases, so also does the width of the plateau. Hence, the time 
necessary for the cage relaxation increases continuously with 
However, in the case of t„ — 3.9 x 10"*, one can observe 
the very beginning of the cage relaxation around i w 2 x 10"*. 
As will be discussed in section|III| this has an important con- 
sequence for the shear rate dependence of ct'"^'''', the stress at 
the maximum of stress-strain curves. 



III. RESULTS AT IMPOSED SHEAR RATE 

An overall shear rate is imposed by moving in the x- 
direction, say, the left wall (zwaii — —20) with a constant ve- 
locity of C/waii- This defines the total shear rate 7tot = CAvaii/^z- 
The motion of the wall is realized in two different ways. One 
method used in our simulations is to move all wall atoms with 
strictly the same velocity. In this case, wall atoms do not 
have any thermal motion. As a consequence, the only way 
to keep the system temperature constant, is to thermostat the 
fluid atoms directly. A different kind of wall motion is realized 
by coupling each wall atom to its equilibrium lattice position 
via a harmonic spring |29]. In this case, the lattice sites are 
moved with a strictly constant velocity while each wall atom 
is allowed to move according to the forces acting upon it [the 
harmonic forces ensure that the wall atoms follow the motion 
of the equilibrium lattice sites]. In such a situation, we can 
thermostat the wall atoms while leaving the fluid particles un- 
perturbed. The temperature of the inner part of the system is 
then a result of the heat exchange with the walls (which now 
act as a heat bath). This method has the advantage of leaving 
the fluid dynamics unperturbed by the thermostat. 

The drawback of thermostating the system through the heat 
exchange with the walls is that, depending on the shear rate 
and the stiffness of the harmonic spring, measured by the 
spring constant fch, a temperature profile can develop across 




FIG. 5: The mean squared displacement (MSD) versus time at tem- 
peratures r = 0.45 (lines) and T = 0.2 (symbols) for various waiting 
times, f„ (increasing from top to bottom), = corresponds to the 
time of temperature change from T = 0.5 to the investigated tem- 
perature. While at T = 0.45 the time translation invariance (TTI) is 
reached after a few thousand LI time units, the data corresponding to 
r = 0.2 indicate an endless evolution towards slower dynamics. The 
larger the wider the plateau and thus the longer the life time of 
the cage. The straight lines are fits to the short time behavior of the 
curves assuming free particle motion with thermal velocity. The ver- 
tical dotted line marks Tco = 2 x 10*. This time is closely related to a 
change in the 7tot-dependence of o-'"^'* [see the discussion of Fig. [SI. 



the system. Note that the smaller the harmonic spring con- 
stant, the better the heat exchange with the walls and thus the 
more efficient the system is thermostated (the imposed shear 
rate having the opposite effect). On the other hand, if fch is too 
small, the fluid particles may penetrate the walls. We find that 
/ch — 25 is a reasonable choice for our model. However, even 
with this value of the harmonic spring constant, we observe a 
temperature profile as the shear rate exceeds 7tot = 10~^. For 
7tot = 10^"^, for example, the maximum temperature in the 
fluid is by about 3% higher than the prescribed value. 

In order to prevent such uncontrolled temperature increases, 
we have therefore decided to apply direct thermostating to the 
inner particles at all shear rates, independently of the possibil- 
ity of the heat exchange with the walls. For this purpose, we 
divide the system into parallel layers of thickness dz — 0.25 
and rescale (once every 10 integration steps) the T/-component 
of the particle velocities within the layer, so as to impose the 
desired temperature T. Such a local treatment is necessary to 
keep a homogeneous temperature profile when flow profiles 
are heterogeneous. To check for a possible influence of the 
thermostat, we compared, for low shear rates (7tot < 10~^), 
these results with the output of a simulation where the inner 
part of the system was unperturbed and the walls were ther- 
mostatted instead. Both methods give identical results, indi- 
cating that the system properties are not affected by the ther- 
mostat. 

However, for wall velocities close to 1 or larger (corre- 
sponding to overall shear rates of 7tot > 2.5 x 10^^), a non- 
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uniform temperature profile develops across the system even 
if the velocities are rescaled extremely frequently |30]. This 
can be rationalized as follows. The heat created by the shear 
motion needs approximately tc = c/Lz to transverse the sys- 
tem (c is the sound velocity). We can estimate the sound ve- 
locity from a knowledge of the shear modulus, G, and the den- 
sity of the system, c = \fGj~P- At T = . 2 we find G « 1 5 (see 
Fig-EJ thus obtaining c« 3.54. A time of tr_ ~ 11.3 is there- 
fore needed for a signal to transverse the whole system. Note 
that the heat creation rate is given by dQ jdt — (T7tot (neglect- 
ing inhomogeneities in the local shear rate). An amount of en- 
ergy equal to ke,T is thus generated within tq — k^T / Q. The 
requirement tq > tc now means that the heat creation must be 
slow enough so that the created energy can be dissipated in the 
whole system efficiently. This gives 7tot < k^T/atc, which, 
after setting and cr«0.6, yields 7tot < 3 x 10"^. 

Figure |6l shows a typical set of (transient) stress-strain 
curves at a temperature of T = 0.2 and for a waiting time 
of = 4 X 10'' LJ time units. The varying parameter is the 
overall shear rate 7tot — U^iaU /Lz (the strain is simply com- 
puted as 7 = 7tott)- First, an elastic regime is observed at small 
shear deformations (7 < 0.02). The stress then increases up 
to a maximum, a^^'^^, before decreasing towards the steady 
state stress at large deformations. Therefore, this maximum 
is sometimes referred to as the yield point |31] or dynamical 
yield stress 1321] . In the following, we will simply refer to this 
quantity as cr^'^'^^, since plastic (irreversible) deformation ac- 
tually sets in before the corresponding value of the strain is 
reached. Moreover, as will be seen below, a'^^'^^ depends on 
strain rate and waiting time in a nontrivial way, so that it is 
difficult, in our simulations, to define a yield stress value from 
such dynamical stress/strain curve. 

As commonly observed in experiments on polymers IjSJ and 
on metallic glasses |^|^, the stress overshoot a^'^^ decreases 
and is observed at smaller strains as the shear rate is lowered 
[see also Fig.lS). Note also that all curves in Fig. |6l show the 
same elastic response at small strains. As also shown in the 
figure, a linear fit to cr = G7 with a shear modulus of G « 15 
describes well the data at small deformations. 

In order to understand the rather strong deviation from lin- 
earity at small strains in the case of 7tot = 10^^, we recall that, 
once the (left) wall starts its motion, a time of approximately 
tc = 11.3 must elapse before the deformation field comprises 
the whole system. This is nicely borne out in the inset of Fig.0 
where, for a wall velocity of U^aii = 0.1, "snap shots" of the 
layer resolved displacement of center of mass (normalized to 
the displacement of the wall) are shown for t — 5 and 11. 
Indeed, the boundary of the deformed region reaches the im- 
mobile wall only after t = 11 LJ time units. We have verified 
this behavior for other wall velocities and have found t « 11 
in all cases. However, as shown in the main part of Fig.0 at 
higher wall velocities, the deformation field is no longer linear 
at the time it reaches the immobile wall. This can be rational- 
ized as follows. The total strain at i = tc is given by 7 = 7tottc 
yielding 7= 11% for 7,01 = 0.4/40 = 10~^. Hence, the elastic 
regime is left already before the whole system is affected by 
the motion of the wall. Putting it the other way, one can esti- 
mate the time for which a locally elastic response can still be 
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FIG. 6: Stress response versus applied strain, 7 = 7iott. for the strain 
rates of 7,0, = 10"^, 10"*, 10"^ and 10"^ Note that both the 
maximum and the steady state values of the stress decrease with de- 
creasing shear rate. All stress-strain curves coincide with the straight 
line (T = G7 (G ~ 15 is the elastic shear modulus) at small deforma- 
tions (7 < 2%). 



observed at a given wall velocity: tei.resp. =7ei/7tot- Assuming 
an elastic response at a strain of a few percent one obtains for 
tei.resp. a time of a few Lennard- Jones units at 7tot = 10^^ [see 
the stars in the inset of Fig.Q. 

The dependence of a^'^'^^ on 7tot is depicted in Fig.|8]for tem- 
peratures of T = 0.2 and T = 0.4. For the lower temperature, 
data are shown for two system sizes — Ly — 10, Lz — ^O 
(averaged over 10 independent runs) and Lx = Ly~ Lz = 40 
(a sole run). As seen from Fig.|SJ for both system sizes, results 
on (jP'^^'' are practically identical. Note that the computation of 
fjpeak ^j. ^^^^ _ 2.5 X 10~^ for the large system required about 
25 days of simulation on a 1.8GHz AMD- Athlon CPU. The 
data point corresponding to 7tot = lO^'' has therefore been 
computed using the average over many small systems only. 
As the results are not sensitive to the system size, we have 
used the smaller system size also in the case of T = 0.4 (again 
averaging over 10 independent runs). 

For T = 0.2, a change in the slope of (TP'^^''-7tot-curve is ob- 
served at a shear rate of approximately 7co = 2.5 x 10~^. 
At shear rates smaller than 7co, the system seems to have 
enough time for a partial release of the stress through rear- 
rangements of particles. Note that the stress overshoot ctP'^^'' 
is observed at strains smaller than 5%. Therefore, small rear- 
rangements are sufficient in order to release the stress consid- 
erably. Indeed, an investigation of the mean squared displace- 
ment shown in Fig. |5lreveals that the MSD departs from the 
plateau for Tco = 2 x 10''. This time is of the same order as 
the inverse of the cross over shear rate thus suggesting that the 
cross over in the 7tot-dependence of o-p'^'''' is related to the be- 
ginning of the cage relaxation. While at higher overall shear 
rates the response of the system is dominated by the (shorter) 
time scale imposed by the shear motion, it is no longer the 
case at 7tot < 7co, where the inherent system dynamics come 
into play. Although not so pronounced, a similar cross over 
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FIG. 7: Short time behavior of the layer resolved displacements 
of the center of mass normalized to the displacement of the moving 
wall, {Xi:m{z; t) — Xcm(z; 0))/ iU^aiii) (z denotes the position of the 
middle of the layer. The system is divided into layers of thickness 
= 1 and Xcm is measured by averaging over the ^-coordinates 
of all particles within the specified layer). The displacement field is 
shown at f = 11 (note that the time needed by the sound to travel 
across the system is given by tc — L^/c ~ 11.3) for various wall 
velocities as indicated in the figure. For Uwaii < 0.2, a linear de- 
formation profile is observed, whereas at higher wall velocities this 
is no longer the case. The inset shows, for a (low) wall velocity of 
U„aii = 0.1, how the deformation field propagates towards the im- 
mobile wall (placed al z — 20). The speed with which the boundary 
of the deformed region extends towards the immobile wall is found 
to be indeed very close to the estimated value of the sound velocity 
c « 3.54. The stars in the inset correspond to U^aii = 0.4 alt = 5 
demonstrating that, at a time corresponding to a smaller strain, the 
local response of the system is elastic [see also the text for more 
discussion]. 

is seen also in the case of T = 0.4 at a larger shear rate in 
agreement with the observation that, compared to T — 0.2, 
the MSD at T = 0.4 leaves the plateau at a shorter time [see 
the MSD(T ^ 0.4) in Fig. [181. Note that, as the structural re- 
laxation time is approximately proportional to the age of the 
system fisll . Tco is of the order of iw The system response 
below the crossover is in fact a complex combination of aging 
dynamics and stress induced relaxation. The aging dynam- 
ics tends to make the system stiffer (see below), so that the 
observed ctP'^^'' is higher than the value one would extrapolate 
from high shear rates. 

The dependence of the stress overshoot a^^^^^ on the im- 
posed shear rate is often expressed with a simple formula 
which goes back to the Ree-Eyring's viscosity theory IT^l34ll . 

a^-'' = <7o + ksT/v*Hj,,,/iyo). (1) 

Here, the activation volume, v*, is interpreted as the charac- 
teristic volume of a region involved in an elementary shear 
motion (hopping) and i/q is the attempt frequency of hopping. 
Obviously, Eq. Q makes sense only at high enough shear 
rates, for in the case of 7tot < i^o, the second term on the right 




I . Q i F 

10"*^ 10"^ 10""^ 10"^ 10 



shear rate 

FIG. 8: The maximum of the stress-strain curve, a^""^, versus strain 
rate. Data for T = 0.2 correspond to two different system sizes: Lx = 
Ly — 10, — 40 (averaged over 10 independent runs) and = 
Ly = Lz— 40 (a sole run). Apparently, results on o-p'^'* do not depend 
much on the system size. At the higher temperature (T = 0.4), only 
the smaller system size (again averaged over 10 independent runs) is 
used. At approximately 7tot ~ 2.5 x 10~^, the slope of cr'"^"'' changes 
significantly at T = 0.2. For T — 0.4, a similar (albeit not so strong) 
change in slope occurs at a higher shear rate. Solid lines are guides 
for the eye. The inset shows the same data, where cr'"^'''' divided by 
temperature is shown versus 7tot. 



hand side of Eq. Q becomes negative. Fitting the data of 
Fig.Hto Eq. Q, we obtain « 2.3 at r = 0.2 and « 3.0 at 
T = 0.4. This result is comparable to the estimates of the free 
volume from experiments on polycarbonate, where a value of 
V* « 3.5nm^ per segment is reported 01. 

Ho Huu and Vu-Khanh have extensively studied the ef- 
fects of physical aging and strain rate on yielding kinetics of 
polycarbonate(PC) for temperatures ranging from — 80°C to 
60°C [note that Tg(PC) « 140°C]. In particular, they have 
measured the tensile stress at yield point, as a function 
of strain rate, e, for various temperatures and different ages 
of the sample. As for the effect of temperature, they find that 
the slope of (T^'(lne)/T (i.e. the activation volume) is prac- 
tically independent of T. Our data also show only a weak 
dependence of v* on temperature, as illustrated in the inset 
of Fig.|8l Note that we have also restricted the data-range to 
higher shear rates where Eq. Q is expected to hold better. 

The above qualitative agreement on the strain rate depen- 
dence of the stress at yield point for our molecular model 
glass and polycarbonate suggests that, for strains smaller than, 
say 10%, the relevant length scale is that of a segment. In 
other words, the chain connectivity has a rather subordinate 
effect on the stress at the yield point (in fact, the connectivity 
becomes important for larger strains, where the well-known 
strain hardening sets in (^|^). 

For the same binary mixture of Lennard- Jones particles as 
in the present work, Rottler and Robbins |33] studied the 
dependence of rj^.^, the maximum of the deviatoric stress, 
on the shear rate. In contrast to our results, no crossover 
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similar to that shown in Fig. |8] was observed in this refer- 
ence. Furthermore, by varying the temperature in the range 
of T G [0.01 0.3] (by a factor of 30), they found that the slope 
of the Tjgy-ln7tot data did practically not change with tem- 
perature, whereas in our case, as discussed above, the slope 
of a^'^'^^-ln 7tot approximately scales with T (see the inset of 
Fig.|8}. Note, however, that in Ref. fsS*] a smaller cutoff radius 
of Tc = 1.5 for the Lennard- Jones potential is used, whereas 
Tc = 2.45 in our model. Furthermore, the pressure in |33] is 
kept at zero at all temperatures, whereas it is always positive 
in our simulations. These differences enhance the repulsive 
(and therefore athermal) character of the system simulated by 
Rottler and Robbins compared to our model. This also ex- 
plains why the shear banding is observed at a temperature as 
low as r = 0.01 in |33], whereas we observe it at T = 0.2 and 
even higher temperatures 12]. It is also worth mentioning that 
the uniaxial strain in Ref. ll33ll was imposed by a simple in- 
stantaneous rescaling of the box dimension and the positions 
of all particles, whereas in our case a more realistic situation 
is considered: The shear strain in the fluid is induced through 
interactions with a moving atomistic wall. We must however 
emphasize that, at the present moment, it is not clear how the 
above differences in details of the model and in the applied 
simulation techniques may lead to the observed discrepancies 
in the behavior of the crP^^'^-ln 7tot curve. 

As an inspection of Fig. |5]reveals, the difference between 
the peak and the steady state stresses decreases as 7tot is re- 
duced thus suggesting that, in the limit of vanishing shear rate, 
^peak converges towards the steady state stress (and therefore 
coincides with the yield stress that could be extracted from ho- 
mogeneous flow experiments). Figure |9] compares these two 
quantities, underlining this expectation further 

It has been shown in experiments on amorphous polymers 
like poly(styrene) and polycarbonate 01|^ that aging strongly 
alters the response of the system to an applied strain. At 
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FIG. 9: The maximum of the stress-strain curves as shown in Fig.|6| 
a^'^, and the steady state stress versus strain rate for a temperature 
of T = 0.2. For a''""^^, we average the results for both system sizes 
shown in Fig.|8| 



small deformations (below, say 5%) the slope of the stress- 
strain curve (elastic shear modulus) increases with progres- 
sive aging. Furthermore, the maximum of the stress-strain 
curve, crP'^^'', is larger for "older" systems and the subsequent 
decrease of the stress ("strain softening" |3]) is more pro- 
nounced. Similar observations are also made in experiments 
on metallic glasses f7]. Interestingly, Fig.llOlshows that these 
features are not limited to polymers or metallic glasses but can 
also occur in simpler models. In Fig.^Jthe stress is depicted 
versus applied strain (defined as 7 = tji^t — tU^aii/ Lz)- Be- 
fore shearing, the system is first equilibrated at a temperature 
of r = 0.5. The motion of the (left) wall is then started at a 
time t„ after the temperature quench. Varying t^, we observe 
similar effects on the stress response as described above. It 
is also observed that, whereas the maximum stress o-p'='* in- 
creases with tv/, the elastic shear modulus (slope of the stress- 
strain curve) seems to saturate already for t„ > 2000 (this is, 
however, hardly distinguishable in the scale of the figure). 

On the other hand, at large deformations, the stress response 
does not show any systematic dependence on the age of the 
system thus indicating a recovery of the time translation in- 
variance: steady shear "stops aging" [17]. In fact, it is well 
known that the shear motion promotes structural relaxation 
and sets an upper bound (^ l/7tot) to the corresponding time 
scale. Once the steady shear state is reached (which is the case 
at deformations comparable to unity), no dependence on the 
system age is expected. Results shown in Fig. ^| are also in 
qualitative agreement with data reported in Ref. |31], where 
the system response to a homogeneous shear was studied via 
Monte Carlo simulations of a binary Lennard-Jones mixture 
(very close to the present model). Note that, in Ref. 1 3 1 ], only 
the contribution to the system response of the so called inher- 
ent structure (configurations corresponding to the minima of 
the energy landscape) has been considered and the effect of 
aging is investigated by applying different cooling rates (not 
by "quenching and waiting" as is the case in our work). De- 
spite these differences in details, results reported in Ref. i3lll 
and our observations are quite similar. More quantitative data 
on the effect of physical aging on the stress at the yield point 
is shown in the inset of Fig. ^| Here, a'?^'^^ is depicted as 
a function of the waiting time, where <„ is varied by more 
than four decades. A logarithmic dependence of a^^'^ on 
is clearly seen for waiting times larger than a few hundred LJ 
time units thus covering about three decades in Such an 
increase in a'^^^^ is consistent with the qualitative idea that the 
system visits deeper energy minima as aging time increases. 
A stronger stress is therefore necessary to overcome the en- 
ergy barriers towards steady flow. It is interesting to note that 
such a ty^ dependence of the stress overshoot is also observed 
in the SGR model |9]. 

As indicated above, simultaneous consideration of figures|8l 
and^Jindicates a rather complex behaviour of o-p'^^'' as a func- 
tion of and 7tot. Considering the similarity in dependence 
for large 7tot or large t^, it is tempting to suggest a rewrit- 
ing of equation^in the form o'P'^^'' = 0-0 + k-^T/ v* ln(7totiw)- 
This modified version of Eq.^does, however, not describe our 
data consistently. At T = 0.2, for example, the (tP'^^''/T versus 
ln(7tottw) curve exhibits different slopes for the data obtained 
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by varying the imposed shear rate (Fig.|8} as compared to the 
simulation results where is the adjustable parameter (cor- 
responding to the data shown in the inset of Fig.llQ>. 

As for the effect of the temperature on the (transient) stress 
response, it is generally known that, due to faster structural re- 
laxation at higher T, the shear stress decreases at higher tem- 
peratures. This is verified in Fig.^Jwhere stress-strain curves 
are shown at T = 0.2, 0.4, 0.43 and 0.5 for a strain rate of 
7tot = 10^'^. Similar to the effect of a decreasing shear rate, 
both the maximum and the steady state values of the stress de- 
crease with increasing temperature. Furthermore, the slope of 
stress-strain curves decreases (the system structure "softens") 
at higher T. Qualitatively similar observations are also made 
on experimental systems (see, for example, figure 1 .20 in 1 12], 
or Refs. OiBElSl)- It is also seen from Fig.^Jthat a change 
of temperature by a factor of two in the glassy state (from 
T = 0.2 to T = 0.4) has less impact on the maximum stress, 
^peak^ than a smaller T-variation close to Tc (from T — 0.4 to 
T = 0.43). This illustrates the sensitivity of the yield point 
to a temperature change in the vicinity of T^. Already from 
this observation, we can expect a similar impact on the T- 
dependence of the static yield stress (see below) close to the 




Strain 

FIG. 10: Aging effects on the stress response to an applied strain. 
(7(7) is shown for a strain rate of 7tot ~ 10~^ at f„ = 20, 2 x 
10^, 4.18 X 10" and 3.99 x lO'^ LJ time units. At small deforma- 
tions, the stress increases faster with progressive aging. The maxi- 
mum observed stress, a""^*, is reached at smaller strains and is higher 
for "older" samples. At large deformations, however, the stress con- 
verges towards the same average value regardless of the age of he 
system. This is a signature of the recovery of time translation in- 
variance, or, equivalently, the erasure of the memory effects due to 
shear induced structural relaxation. The inset shows the variation 
of the maximum stress with the waiting time. Note that, here, tw 
is varied by more than 4 orders of magnitude, i.e from t„ = 20 to 
tw ~ 3.99 X 10"". The solid line is a guide for the eye. The data cor- 
respond to a system size of = Ly = — 40 (a sole run). For the 
largest waiting time [t„ = 3.99x 10^ (=2x lO^MD steps)], however, 
average over 10 independent runs of a smaller system size is used 
{Lx = Ly — 10, Lz — 40). Note that already at this (smaller) system 
size, the size effects are practically negligible [see also Fig.|8|. 



mode coupling critical temperature of the system. 



IV. RESULTS AT IMPOSED STRESS 

In this section we study the response of the system to im- 
posed shear stress. The system is prepared in a similar way 
as described in previous sections so that, at the beginning of 
the measurement, the structural relaxation times of the system 
are much larger than the time scale of the simulation. Starting 
with a{t = 0) = 0, we gradually increase the external stress 
(i.e. the force acting on the atoms of the left wall) and record 
quantities of interest, such as the internal energy, the stress 
across the system, the center of mass velocity of the walls and 
of the fluid, etc... 

It is generally accepted that imposing an external stress 
leads to a shift in the density of accessible states towards 
higher energy configurations. For the binary Lennard-Jones 
model of the present work, Fig.[21shows the potential energy 
per particle, Cpot, as measured in simulations where the im- 
posed stress is periodically varied the range a G [—0.8 0.8] 
(see the zigzag line in Fig. [12] Similar stress ramps were also 
used by He and Robbins |29] in order to determine the static 
friction between two solid bodies mediated by a layer of ad- 




strain 

FIG. II: Effect of temperature on the stress response to an ap- 
plied strain. cr(7) is shown for a strain rate of 7tot = 10""^ at 
T = 0.2, 0.4, 0.43 and 0.5 (note that the mode coupling critical 
temperature of the system is Tc — 0.435). For temperatures below Tc 
the system was first aged during t„ = 4 x 10* LJ time units before the 
beginning of the measurement. Similar to the effect of a decreasing 
shear rate, both the maximum and the steady state values of the stress 
decrease with increasing temperature. Furthermore, at small strains, 
the slope of stress-strain curves (elastic shear modulus) decreases in- 
creasing temperature (see the inset) thus indicating a softening of the 
system structure at higher T [compare to Fig.|5|. Note, however, that 
these changes are much more pronounced in a narrow temperature 
interval around Tc. Results here correspond to averages over 10 in- 
dependent runs. The system size was = Ly = 10 and Lz ~ 40. 
The inset shows a magnification of the small strain region of the same 
data. 
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sorbed molecules). 

Note that the maxima and minima of the potential energy 
correspond to \<j\ = 0.8 and a — Q respectively. Starting at a 
minimum of epot (cr = 0), the potential energy fluctuates for 
a while around this minimum before increasing sharply to- 
wards a maximal value. This corresponds to a branch where 
\<j\ increases from to 0.8. The descent from this maximum 
towards the subsequent minimum (|cr| decreases from 0.8 to 
0) is, however, more gradual and indicates a dependence of 
Cpot on the stress history. Finally, we also observe that, at 
high (7, the quiescent energy distribution observed at small 
stresses at the very beginning of the stress ramp simulation, 
is never reached again whereas the stress itself passes through 
zero periodically. This dependence on &, however, is con- 
siderably weakened as the stress increase rate reaches values 
below 5 X 10^^. 

While the potential energy per particle is easily measured 
in a simulation, this is not the case in real experiments. The 
velocity of the solid boundary (upon which the stress acts), 
however, is experimentally accessible. Fi gure [T3l depicts the 
wall velocity measured in simulations at = 5 x 10"^. Fol- 
lowing the convention, the applied stress is shown on the ver- 
tical axis, whereas on the horizontal axis the system response 
is depicted. We first note that, at small stresses, the system 
resists to the imposed stress and thus prevents the wall from 
moving. Only when the magnitude of the stress exceeds a cer- 
tain (yield) value, a non-vanishing wall velocity is observed. 
Furthermore, after a cross over regime around the threshold 
value of the stress, the wall velocity increases almost linearly 
with stress increment. 
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FIG. 12: Effect of the rate of stress increase on potential energy 
per particle, e^at is measured during cyclic variations of the imposed 
stress as sketched by the zigzag hue (note that a varies in the range 
[—0.8 0.8]). The horizontal axis counts the number of cycles. Three 
rates of stress variation are shown here, a = 5 x 10^^, 5 x 10~^ 
and 10~^ . T he higher a, the higher the potential energy per particle 
at small stresses (minimum of epot). The vertical dashed lines mark 
simultaneously the maxima of epot and Icrj. They serve to better rec- 
ognize the asymmetry of epot on both sides of the stress maximum 
and recall the presence of a hysteresis effect. 



On the other hand, as the magnitude of the stress is de- 
creased again, the wall motion first slows down along the same 
line as in the stress increase case but then departs towards 
higher wall velocities. A hysteresis loop is thus formed as 
expected from an analysis of the asymmetry of epot around the 
stress maximum [see Fig. 1121 . Similar observations are made 
in experiments on pastes, glass beads, dense colloidal suspen- 
sions 1 7] and foams |8]. Note also that, as expected from the 
symmetry of the system response with respect to positive and 
negative stresses, the shape of the observed hysteresis loop is 
identical for both directions (signs) of the applied stress. 

Next, we investigate the dependence of the system response 
to an applied stress on the rate of stress variation. For this pur- 
pose, a is varied by two orders of magnitude, from 5 x 10^^ 
to 5 X 10~^. Figure IT?! depicts stress ramp data now aver- 
aged using the symmetry with respect to negative and positive 
stresses. Again, for all values of & shown in this figure, no 
flow is observed for too small stresses (below, say 0.4). How- 
ever, for a given stress above, say, a = 0.7, the wall velocity 
is lower at higher &. To put it the other way, when |cr| is in- 
creased faster, a given wall velocity is reached at a higher \a\, 
i.e. on a later time. This may be rationalized by noting that, 
at a higher stress increase rate, the system has less time to de- 
velop a response corresponding to the actual (instantaneous) 
stress. Therefore, the mobility increase corresponding to an 
increase of the stress is retarded and is observed later, i.e. at 
higher stress. 

However, it is also seen from Fig. O that, akeady at fi < 
2 X 10^^, the effect of a on the system response is of order 
of the measurement uncertainty, so that no systematic depen- 
dence on a can be seen for tr < 2 x 10~^. This is consistent 
with the behavior of the potential energy per particle which 
becomes practically independent of & in the same (j-range [see 
Fig. 1121 . Therefore, we may describe this regime of slow vari- 
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FIG. 13: The applied shear stress (vertical axis) and the resulting 
wall velocity (horizontal axis) measured during stress ramps with a 
rate of a = 5 x 10" ''. The result shown here is an average over two 
independent runs each containing 15 full cycles of stress variation 
[see the zigzag line in Fig. ll2l . 
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ation of cr as quasistatic. 

Results presented above and in previous works 
show that our model system shares many features of the so 
called soft glassy materials. In particular, the existence of a 
yield stress is suggested in Figs. I13landll4l Figure dis- 
plays further evidence of the existence of a yield stress still as 
the dramatic change in the wall velocity at a threshold stress 
value is emphasized using a logarithmic scale for the horizon- 
tal axis. In a narrow stress range around a ~ 0.6, the wall ve- 
locity and thus the overall shear rate increases approximately 
by three orders of magnitude [see also Fig.^. Again, a linear 
regime is observed at high stresses f?"]. Besides the hystere- 
sis already discussed above, an investigation of the decreasing 
branch on the stress-wall velocity curve in Fig.llSlreveals that, 
as (7 falls below a certain value, the wall velocity becomes 
even negative [see the inset]. This clearly illustrates the pres- 
ence of attractive forces which, now, are stronger than the im- 
posed stress and thus capable of reducing the amount of strain. 
Indeed, an inspection of the center of mass position of the wall 
and of the fluid shows that both these quantities exhibit a max- 
imum at the place where the velocity passed through zero (as 
the stress is further reduced, X^^ decreases in accordance with 
the observation of a negative velocity). Very similar observa- 
tions are also reported on the experimental side ^TJ . 




wall velocity 



FIG. 15: The applied shear stress (vertical axis) and the resulting 
wall velocity (horizontal axis) measured during stress ramps with a 
rate of a = 2 x 10^'' . The vertical solid line roughly marks the elastic 
contribution to the strain rate (7^' ~cr/G. Elastic deformation gives 
rise to a non vanishing wall velocity even at the smallest imposed 
stress. See also Fig.|S|for an estimation of the shear modulus G). The 
inset shows that, as the stress is decreased, the wall velocity changes 
its sign thus indicating that the attractive forces are stronger than 
the imposed stress so that the direction of deformation is reversed in 
order to reduce the amount of the accumulated strain. 
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FIG. 14: The hysteresis loops as measured during stress ramps 
with rates of stress variation of a = 5 x 10~® (filled diamonds). 
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The (T = 5 X 10~ -curve is an average over 10 independent runs with 
a unique variation of the stress from to 0.76. The remaining curves 
correspond to averages over two independent runs each containing 
many full cycles of stress variation in the intervall [—0.8 0.8]. The 
innermost loop (filled diamonds) corresponds to the smallest cr. Note 
that the surface of the hysteresis loop increases at higher stress varia- 
tion rates thus indicating stronger retardation effects. Note also that, 
for the two highest ct, the loop does not close within the simulated 
stress range. It would close at much higher stresses than shown in 
the figure. 



V. MEASUREMENT OF THE YIELD STRESS 

As discussed in section|IJ there seems to be a close connec- 
tion between the existence of a yield stress and the observation 
of the shear banding phenomenon in many soft glassy materi- 
als 0,0 HH. In particular, it is commonly expected that, in a 
state where the yield stress vanishes (at high temperatures, for 
example) the shear bands should also disappear, i.e. the whole 
system should flow. In addition to this experimental aspect, a 
study of the yield stress is also motivated from the theoretical 
point of view. For example, the so called soft glassy rheol- 
ogy model (SGR) of Sollich issll (an extension of the trap 
model f3^ taking into account yielding effects due to an ex- 
ternal flow) predicts a linear onset of the dynamic yield stress 
as the glass transition is approached, cr(7tot — > 0) ^ 1 — x. 
Here, x is a noise temperature, x= \ corresponds to the glass 
transition (or "jamming") temperature, and x < \ character- 
izes the glassy or "jammed" phase. On the other hand, nu- 
merical studies of a p-spin mean field Hamiltonian 1 16] pre- 
dict that the dynamic yield stress vanishes at all temperatures. 
There has recently been a more microscopic approach based 
on an extension to non equilibrium situation 1 37] of the mode 
coupling theory of the glass transition (MCT) |38]. An anal- 
ysis of schematic models within this approach shows a rather 
discontinuous change in the dynamic yield stress at the mode 
coupling critical temperature, T^. 

The reader may have noticed that the above mentioned the- 
ories make predictions on the dynamic yield stress [defined 
as (7(7101 0)]. Our interpretation of the shear banding, 
however, makes use of the idea of resistance to an applied 
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stress which is related to the presence of a static yield stress. 
Similar to the difference between the dynamic and static fric- 
tion isoll . the static and the dynamic yield stresses are not nec- 
essarily identical. Indeed, for our model glass, we find that 
CTy > cr(7tot — > 0) [see Fig.lJJ. Therefore, a measurement 
of the static yield stress gives at least an upper bound for the 
dynamic counterpart. As we will see below, the static yield 
stress decreases rather sharply as the mode coupling critical 
temperature of the model {T^ — 0.435) is approached. Unfor- 
tunately, when measuring Cy at temperatures close to T^, one 
is faced with the problem that the time scale imposed by the 
external force (which if of order of the inverse stress variation 
rate, i.e. — cr/a) and that of the (inherent) structural relax- 
ation, Treiax, are not well separated. In particular, the condition 
ta <C Treiax is not Valid at temperatures close to Tc. Therefore, 
as will be discussed below in more details, a conclusive state- 
ment on the interesting limit of (Jy{T T^) can still not be 
made. 

Preliminary results on the static yield stress have been re- 
cently obtained within the driven mean field p-spin mod- 
els |4ul . Using the fact that the free energy barriers are finite at 
finite system size, the model has been investigated by Monte 
Carlo simulations in the case of p — ?> for a finite number 
of spins, thus allowing the thermal activations to play a role 
which they could not play in the case of an infinite system size. 
Results of these simulations support the existence of a critical 
driving force below which the system is trapped ('solid') and 
above which it is not ('liquid') |40]. Results based on this new 
approach on the temperature dependence of the yield stress 
and, in particular, on its behavior close to Tc are, however, 
lacking at the moment. 

Here, we adopt a method very close to a determination of 
the (static) yield stress in experiments, i.e. we use the defini- 
tion of CTy as the smallest stress at which a flow in the system 
is observed. As we are interested in a study of the tempera- 
ture dependence of CTy and, in particular, in CTy (T) close to the 
mode coupling critical temperature, we have varied the tem- 
perature in the range of T G [0.1 0.44] (recaU that Tc = 0.435). 
For each temperature, ct was increased stepwise by an amount 
of da — 0.02 once in each 1000 LJ time units during which the 
velocity profile corresponding to the imposed stress is mea- 
sured. Among other quantities, we also monitor the motion 
of the center of mass of the wall and also of the fluid itself. 
Note that the overall stress increase rate in these simulations 
is CT = 2 X 10^^, and thus corresponds to a quasi static varia- 
tion of the stress [see the discussion of Figs.[T2land [T4l . For 
each temperature, the simulation was performed using 10 in- 
dependent initial configurations. 

Recall that there is always an elastic contribution to the sys- 
tem response to an applied stress. The corresponding center 
of mass velocity can simply be estimated as V^^ — a/G. This 
contribution is negligible at lower T for two reasons: (i) due 
to the high stiffness of the system (large G), V^}^ is relatively 
small and (ii) the onset of the shear motion is quite sharp at 
low T thus leading to much higher velocities (compared to 
as soon as the applied stress exceeds CTy. In contrast, 
close to Tc, the shear modulus is quite small [see, for exam- 
ple, the slope of the stress-strain curve at T = 0.43 in Fig. II II 



thus leading to a larger 1/,'^. Furthermore, there is no sharp 
variation in V^m as a function of applied stress. For a mea- 
surement of CTy close to Tc, it is therefore important to correct 
for the elastic contribution to the system response. For this 
purpose, we have determined the T-dependence of the shear 
modulus. The center of mass velocity of the fluid has then 
been corrected subtracting, for each temperature, the corre- 
sponding V^^ = ct/G. 

Figure [^depicts the applied stress (vertical axis) and the 
resulting (corrected) center of mass velocity of the fluid, X^m, 
averaged over all independent runs (horizontal axis). A log- 
log plot is used in order to emphasize the continuous variation 
of Vcm with decreasing stress at high temperatures. Contrary 
to low temperatures (T < 0.35) where a plateau followed by a 
sharp drop towards zero in V^m is observed, the center of mass 
velocity of the fluid at high temperatures decreases rather con- 
tinuously for small stresses. 

As a first attempt to determine the yield stress, we apply 
linear fits to the data shown in Fig.^] As shown in the same 
figure, the chosen fit range roughly corresponds to the plateau 
region at low temperatures. For T < 0.35, we thus expect the 
fit result not to be significantly different from the "real" value 
of CTy. However, as an investigation of the high-T behavior of 
Vcm in Fig. ^] suggests, this method is not expected to give 
accurate results for CTy at high temperatures (T > 0.38, say). 

A slightly different approach in determining CTy is to find the 
smallest stress for which the center of mass velocity exceeds 
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FIG. 16: Effect of the temperature on the response of the system to 
imposed shear stress. The imposed stress is shown on the vertical 
axis while the center of mass velocity of the fluid (inner part of the 
system) is depicted as horizontal axis. Each curve corresponds to an 
average over 10 independent runs. The stress was increased stepwise 
by an amount of da — 0.02 once in each dt = 1000 LJ time units 
((T = 2 X lO^''). Note that the contribution of the elastic deformation 
to the center of mass velocity (V^^ ~ Lz&/{2G), using the value of 
G~15 atT = 0.2) has already been subtracted from the data. Note 
also that the statistical uncertainty of Vcm is approximately of order 
of 10~*. The vertical dashed lines show the limits of the Km-range 
used in the fit to cr = cry + aVcm- The inset is a magnification of the 
high stress regime. 
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a certain, small value, V^™". Here, we further require that V^m 
must remain larger than V™^ for all subsequent stresses. This 
last condition serves to reduce errors due to fluctuations of 
VI;tn. In applying this definition, we use the result of each inde- 
pendent run on Va-n separately and thus obtain, for each VJ,™", 
a set of yield stress values. This allows an estimate of the sta- 
tistical error. Figure^]compares the yield stress obtained via 
the linear fit to Vcm with results of the second approach for 
Km" = lO^^'i 10"^ ™d 10"^. Not unexpectedly, it is seen 
from Fig.^]that the quality of results on cTy strongly depends 
on temperature. At temperature far enough from Tc, say, for 
T < 0.35, CTy is rather insensitive to a change of V^™" (dif- 
ferences caused by various choices of V^™" are of order of the 
statistical error). At higher temperatures, however, the varia- 
tion of CTy with the choice of V^™" is remarquable: At T = 0.43 
it varies between 0.14, 0.22 and 0.33 for y™" = 10"*, 10"^ 
and 10~^. Therefore, result on CTy at temperatures close to Tc 
should be considered as rough estimates only. 

The origin of the difficulty in estimating the static yield 
stress of the system at temperatures close to Tc, can be un- 
derstood by comparing the time scales relevant to the prob- 
lem. First, there is a time scale related to the imposed stress 
tcr = ct/ct. The second relevant time scale is that of the struc- 
tural relaxation, Tieiax. The static yield stress is well defined 
in the limit of a quasi static variation of stress, i.e. ct 
(to- oo) while at the same time keeping r^iax ^ t&- Us- 
ing CT « 0.5 and ct = 2 x 10^^ (note that this value of ct was 
used at all temperatures in order to determine the yield stress) 
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FIG. 17: The effect of the temperature on the static yield stress, 
(Ty. The solid line shows cry obtained from fits to ct = o-y + aVcm {a 
and CTy are fit parameters) using the fit range Vcm G [0.0003 0.005] 
[see Fig. 1161 . The symbols correspond to CTy defined as the small- 
est stress for which (and for all subsequent higher stresses) the wall 
velocity exceeds a certain minimum value, V^™°. Three choices of 
V"cm" are compared: 10"" (circles), 10"^ (diamonds) and 10"^ (tri- 
angles). While CTy is relatively insensitive to a choice of V^™" at low 
temperatures, it is not the case for temperatures close to Tc, where 
it continuously decreases as the V^"'" is reduced. The vertical arrow 
marks the mode coupling critical temperature Tc = 0.435. For clarity, 
error bars are shown for the case of V^™" — 10"'' only. 



we obtain ~ 2 x 10^. We are therefore led to verify if the 
condition T^eiax ^ 2 x 10** is satisfied at all temperatures. For 
this purpose, we define Treiax as the time needed by the mean 
square displacement of a tagged particle to reach the particle 
size. Figure ^] shows the mean square displacement of the 
unsheared system for T e [0.1 0.44] (recall that = 0.435). 
For all these temperatures, the waiting time between the tem- 
perature quench (from an initial temperature of T = 0.5 to 
the actual temperature) and the beginning of the measurement 
was t„ = 1-8 X 10"*. At low temperatures, the MSD prac- 
tically remains on a plateau for the whole duration of the 
simulation indicating that Treiax is much larger than the sim- 
ulated time of 2 X lO'* LJ time units. At higher temperatures 
(T > 0.41), however, after a long plateau, the MSD eventu- 
ally enters the diffusive regime and reaches a value compara- 
ble to unity within the simulated time window. Obviously the 
condition Treiax ^ ^<t IS violated at these temperatures. Hence 
at least for a waiting time of t„ = 4 x 10* and for the choice 
ofCT = 2xlO"^, the computed static yield stress is not well 
defined close to Tc. 



VI. CONCLUSION 

Results on the yield behavior of a model glass (a 80:20 bi- 
nary Lennard-Jones mixtures studied by means of molec- 
ular dynamics simulations, have been reported. One of the 
major motivations of the present work is the observation of 
shear localization (below the glass transition temperature and 
at low shear rates) in the present model and the suggestion of 




FIG. 18: The mean square displacement (MSD) of a tagged parti- 
cle averaged over all three spatial directions at various temperatures 
ranging from the supercooled state (T = 0.44 > Tc = 0.435) down 
to the frozen state T = 0.1. The time tw indicates the time elapsed 
between the temperature quench and the beginning of the measure- 
ment. At a time of t = 2 x 10*, the MSD hardly leaves the plateau at 
low temperatures. For temperatures T > 0.41, however, it approx- 
imately reaches the size of a particle within the same time interval 
indicating that a complete structural relaxation has taken place. The 
horizontal dotted line marks MSD = 1. 
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a link between this phenomenon and the existence of a static 
yield stress |2] (under applied stress, the system does not flow 
until the stress exceeds a threshold value). A particular em- 
phasis thus lies on the yield stress and its dependence on tem- 
perature. 

First, the system stress-strain curve under startup of steady 
shear has been studied. The effect of physical aging (charac- 
terized by the waiting time, <„), shear rate (7tot), and temper- 
ature on the stress-strain relation has been investigated. Re- 
gardless of these parameteres, all observed stress-strain curves 
first exhibit an elastic regime at small shear deformations 
(7 < 0.02). The stress then increases up to a maximum, ctP'^*, 
before decreasing towards the steady state stress at large de- 
formations. The steady state stress (corresponding to large 
deformations) shows a dependence on temperature and on the 
applied shear rate, but is independent of the system history, 
indicating a recovery of the time translation invariance due 
to shear induced structural relaxation fW\. In contrast, the 
stress overshoot a'^'^^^, (the first maximum of the stress-strain 
curves, spmetimes described as a dynamical yield stress) de- 
pends on the imposed shear rate and on the waiting time 
(physical aging). It is observed that, at relatively high shear 
rates or for large waiting times, the maximum stress increases 
with ln(7tot) or with \n{t^), respectively. These observations 
are consistent with experiments on amorphous polymers lOlQ] 
and on metallic glasses |5, 6|, and also correspond to the be- 
haviour predicted using the soft glassy rheology model [9] . 

For shear rates below a certain, cross over shear rate, 7co, 
however, a decrease in the slope of cr'"^^''-7tot curve is seen. 
A comparison with the steady state shear stress suggests that 
^peak saturates at the steady state stress level as the imposed 
shear rate approaches zero. Moreover, an analysis of the 
mean square displacements of the unsheared system reveals 
that the cross over shear rate, 7co, is very close to 1/tco, 
where t^o marks the time for which the mean square displace- 
ment gradually departs from the plateau-regime [see Fig.lSj. 
We therefore associate this crossover with the beginning of 
the cage relaxation, which leads to the possibility of small 
(compared to the size of a particle) rearrangements thus al- 
lowing at least a partial release of the stress, gammadotco 
is also comparable to the inverse of the waiting time: for 
gammadottot < gammadotco, the response of the system 
is directly influenced by the aging dynamics. 

In order to build a closer connection between our studies 
and typical rheological experiments, stress ramp simulations 
are performed and the system response is analyzed for stress 
increase rates ranging from a — 5 x 10^^ to <t = 5 x 10^^. In 
agreement with experiments on complex systems like pastes, 
dense colloidal suspensions ITd and foams hysteresis 
loops in the system response are observed. These loops be- 
come wider as a increases. An analysis of the potential en- 
ergy per particle for different a nicely shows how high energy 
configurations are favored by the faster stress variations. This 
also yields an estimate of quasi static stress application. We 



find that, for our model, = 2 x 10~^ is slow enough so that 
simulations with this stress variation rate can be used in order 
to obtain a reliable estimate of the static yield stress. 

Finally, the static yield stress, ay, is determined and its re- 
liability is discussed. Our numerical results confirm the ob- 
servation of reference (2], that the static yield stress is higher 
than the low shear rate limit (7(7 0) observed in steady 
shear experiments. The system can therefore produce shear 
bands for stresses in the range [a{j — > 0),f7y]. 

At temperatures far below the mode coupling critical tem- 
perature of the model (T^ — 0.435), a slight increase of dy 
with further cooling is observed. At temperatures close to T^, 
however, the static yield stress strongly decreases as T is in- 
creased towards T^. As to the reliability of the data, relatively 
accurate estimate of ay is obtained at low temperatures (for 
T < 0.35). Results on the yield stress at temperatures close 
to Tc, however, are very sensitive to the applied criterion. An 
investigation of the dynamics of the unperturbed system re- 
veals that, for T close to T^, the structural relaxation times are 
far from being large compared to the time scale imposed by 
the external force (the inverse of the stress increase rate, a/&). 
Therefore, for the simulated waiting time of 4 x lO"*, the static 
yield stress is no longer well defined at these high tempera- 
tures. This underlines the fact that a very good separation of 
time scales between the experimental and intrinsic time scales 
is necessary in order to properly define a static yield stress. 

It must, however, be emphasized that, even though an in- 
crease of & apparently leads to a validity of Tjeiax ^ t^, this 
would violate the condition of a quasi static variation of the 
stress. A more physical way to improve the accuracy of re- 
sults on ay is to increase the waiting time, in order to allow 
Treiax to grow bcyond t^. Noting that, at higher temperatures 
(but still below Tc), Treiax increases less strongly with t^, (inter- 
rupted aging), the limit of large Treiax becomes progressively 
more time consuming in terms of computation time. 

Our numerical study shows that a very simple model, stud- 
ied numerically on relatively short time scales, can exhibit 
most of the complex rheological behaviour of soft glassy sys- 
tems, but also of "hard" (metallic) glasses (it is interesting in 
this respect to note that the simulated system was originally 
intended to mimic a NiPd metallic glass). This suggests that 
these features are generic to most glassy systems, although in 
practice the values of the parameters may considerably vary 
from system to system. 
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